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Abstract. We study excitation energy transfer in a simple model of photosynthetic 
complex. The model, described by Lindblad equation, consists of pigments interacting 
via dipole-dipole interaction. Overlapping of pigments induces an on-site energy 
disorder, providing a mechanism for blocking the excitation transfer. Based on the 
average efficiency as well as robustness of random configurations of pigments, we 
calculate the optimal number of pigments that should be enclosed in a pigment- 
protein complex of a given size. The results suggest that a large fraction of pigment 
configurations are efficient as well as robust if the number of pigments is properly 
chosen. We compare optimal results of the model to the structure of pigment-protein 
complexes as found in nature, finding good agreement. 



1. Introduction 

Photosynthesis is the main natural process for harvesting Sun's energy on Earth, 
providing a food source for a great variety of organisms ranging from highly evolved 
photosynthetic systems in higher plants to simpler bacteria pQ . In spite of such diversity, 
basic underlying principles are shared among majority of light harvesting organisms. 
The initial stage of photosynthesis usually involves multiple pigment protein complexes 
(PPC) that consist of a number of pigment molecules (i.e., chromophores) held in place 
by a protein cage. PPCs are employed to absorb incoming photons as well as to transport 
the resulting excitation to the reaction center, where the excitation is used to initiate 
chemical reactions. The absorption of light and in particular channeling of the absorbed 
energy to the reaction center is known to achieve high efficiency pQ. 

One of the most studied PPCs is the Fenna-Matthews-Olson (FMO) complex that 
is found in the photosynthetic apparatus of green sulfur bacteria. It is the first PPC 
for which the atomic structure has been determined [3J. FMO is composed of a large 
protein envelope that encloses a tightly packed group of 7 pigment molecules^] called 
bacteriochlorophylls (BChl). In FMO the main role of BChl pigments is actually 
not to absorb light but instead to transport electronic excitations from the input 
pigment, which is close to the antenna complex, towards a "sink" pigment channeling 

| Recent structural analysis [3] suggests also the presence of an eighth BChl that is weakly bound to 
each monomeric unit as an additional input site. Its distance from the core 7 pigments is quite large 
and therefore we do not take it into account in our study. 
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excitation to the reaction center. Recently discovered long-lasting quantum coherence in 
FMO complex [3] gave an additional boost to studies of excitation transport in PPCs. 
Previously it was namely believed that the transport in PPCs is predominantly of a 
classical nature. Many aspects of time dependence have been studied [SI El El, including 
the functional role of quantum coherences [SJ EJ QUI E], as well as the possibility of 
transport enhancement via environmental interaction [TTJ Q21 Q21 03] ■ 

Also important is the question of structural characteristics of PPCs that enable 
efficient excitation transfer, for instance, why has a given complex precisely the shape 
found in nature. Positions and orientations of pigments, known with high precision 
from crystallographic measurements, on casual inspection show no clear ordering or 
organization that would enable an easy classification of efficient configurations. It is 
also not clear, how special efficient configurations are - whether efficient configurations 
are a result of a long-lasting process of evolutionary improvement, or, can efficient 
configurations be readily achieved probing few random configurations of pigments. Prior 
to the availability of crystallographic structural data, the average minimal distance 
between pigment molecules was estimated based on the comparison between florescence 
yield of in vivo and in vitro chlorophyll solutions [32] . Recently, the efficiency of random 
configurations within simplified models of PPC was inspected [TBI H3 UHl HH] , suggesting 
that the efficient configurations are relatively probable. This complies with the results 
obtained for the Photosystem I from plants and cyanobacteria [2TJ1 12"T] , where random 
orientations of pigment molecules were probed, and high efficiencies of excitation energy 
transfer (EET) were obtained irrespectively of the pigment orientation. Also, the PPC 
configuration was shown to be robust to the removal of a pigment from the complex 

Fundamental, yet still unanswered question that we address in present work is, why 
has a particular photosynthetic complex exactly the specified number of pigments. In 
other words, what is the optimal number of pigments for a given size of the complex^ or, 
equivalently, what is the optimal size of the complex for a given number of pigments. For 
instance, why do we find exactly 7 pigments in the FMO complex and not more or less. 
Using a simple model whose parameters are taken from experiments, that is without 
any fitting parameters, we calculate the optimal number of photosynthetic pigment 
molecules for different complex sizes and compare these theoretical predictions with the 
actual number of pigments found in naturally occurring complexes. We find a very good 
agreement for a variety of PPCs in different organisms. To judge the optimality we use 
two criteria: (i) the average efficiency of the excitation transfer, where the averaging is 
performed over random configurations of pigment molecules, and (ii) the robustness of 
the efficiency to small variations of pigment's locations. A rationale behind these two 
choices is that "good" PPCs should have high efficiency but at the same time also be 
robust. A specially "tuned" configuration, that has a very high efficiency which though 
is very fragile, will obviously not work in a natural environment with its changing 



§ Note that the size of the PPC might be fixed by external factors like for instance the membrane 
thickness. 
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conditions. Additionally, from an evolutionary perspective, the efficiency should be 
stable with respect to different foldings of the protein cage. It is advantageous to have 
a PPC with such a number of pigments that will results in high average efficiency, 

1. e., in many close-to-optimal pigment configurations. Thereby, a small change in the 
environment, be it of a chemical origin or for instance a genetic mutation, will still 
result in a functional PPC. Robustness of quantum coherence to structural changes in 
the PPCs has been also found experimentally [22] • The model we use to describe the 
excitation transfer across the PPC consists of a Lindblad master equation describing 
a dipole-coupled pigments with an on-site excitonic energies being determined by the 
distances between disc- like pigments. If two discs come too close, i.e., if they overlap, this 
effectively rescales their energies, introducing a disorder. We should say that the optimal 
number of pigments that we predict is quite insensitive to details of the underlying 
model. Agreement between predictions of our model and naturally-occurring PPCs 
shows that Nature has optimized PPCs by using just the right number of pigments so 
that the resulting PPCs are highly efficient and robust at the same time. 

2. The model 

To calculate the efficiency of the excitation transfer in the PPC, we need equations 
of motion describing dynamics of excitation on multiple chromophores, coupled to the 
environment. It is the ratio of chromophore-chromophore interaction strength to the 
chromophore-environment coupling that determines the applicability of various models. 
When interaction between chromophores is small compared to the environmental 
coupling the Forster theory [23] is applicable, leading to a picture of incoherent hopping 
of excitation between chromophores. In the opposite limit of strong inter-chromophore 
interaction and weak coupling to the environment, the excitation dynamics can be 
described by quantum master equation, either the Redfield equation or, employing 
a secular approximation, the Lindblad equation |24J. For nonperturbative parameter 
ranges, more advanced methods have been developed [2SJ [251 12Z] , usually at the expense 
of higher computational complexity. 

Our main optimality criterion, namely the transfer efficiency, is relatively 
insensitive to details of excitation time-dependence. Thus, we are going to use the 
simplest description of EET with the Lindblad equation, which retains coherent nature 
of transport, while still taking environmental interactions into account. We expect that 
more exact descriptions, which in general enhance excitonic oscillations, lead to similar 
results due to our averaging procedure. Also, these oscillations appear on the time scales 
of few 100 fs, which is much shorter than the time scale of excitonic transfer. 

In the following, we will introduce the Lindblad equation for the overlapping disc 
model, used for the description of excitation dynamics in PPCs. Optimality criteria for 
the efficiency of PPC configurations as used in latter sections will also be presented. 
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2.1. Lindblad master equation 

Internal dynamics of the system of N chromophores within a single-excitation manifold 
is determined by the Hamiltonian of the form [2H] 

N N 

H = J] e n \n)(n\ + V mn(\m)(n\ + \n)(m\), (1) 

n=l n^=m=l 

where a state \n) represents an excitation on the n-th chromophore site, i.e., the 
electronic state of the n-th chromophore being in the 1st excited state. Because EET 
is sufficiently fast, events with two excitations being present at the same time are rare 
and it is sufficient to consider only zero and single-excitation subspace [H] . The coupling 
V mn is due to dipole-dipole interaction between chromophores of the form 

^" '-U \ mn mn 



where r mn vector, connecting the m-th and n-th chromophores, d n is 

a transition dipole moment between the ground and the 1st excited state of the n-th 
chromophore. 

Because the system of chromophores is coupled to the protein and nuclear degrees of 
freedom it is described by a reduced density matrix p. Decoherence due to environmental 
interaction, recombination of excitation to the ground state, and transfer of excitation 
to the sink, are modeled by Lindblad superoperators that augment the von Neumann 
equation for the time evolution of density matrix, 

P = ~i[H, p] + £dcph(p) + £sink(p) + £rccomb(p)- (3) 

To model the effects of the environment, we have taken a simplified picture of purely 
dephasing Lindblad superoperators (i.e. Haken-Strobl model), which is believed to 
capture the basic environmental effects and was used in various previous studies [ZS1 
[TTj [T3"| IT2] . For longer time, relevant for the efficiency of PPC, it has been shown 
that the description with the Lindblad equation accounts for the main features of the 
dynamics [301 EI] • Dephasing Lindblad superoperator destroys a phase coherence of any 
coherent superposition of excitations at different chromophores, and is given by 

N ( 1 \ 

Aieph(p) = 27^ (\n)(n\p\n)(n\ - -{|n)(n|, p} ) , (4) 

n=l V L J 

where a site-independent dephasing rate is given by 7 and { , } represents the 
anticommutator. Irreversible transfer of excitation from the iV-th chromophore to the 
sink \s) is modeled by Lindblad superoperator 

£si„k(p) = 2k (\s)(N\p\N)(s\ - ^{\N)(Nlp}^j , (5) 

where k denotes the sink rate. The irreversible loss of excitation due to recombination 
is given by an analogous term 

- / 1 \ 

£™b(p) = 2T £ |0)(n] p |n)(0| - -{\n)(n\,p} , (6) 

n=l V Z / 
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with a site-independent recombination rate given by Y. The ground state of a 
chromophore system (state without any excitation) is represented as |0). Note that 
£ S i n k and Aecomb can be equivalently represented by an antihermitian Hamiltonian at 
the expense of non-conserved density matrix probability [12], avoiding the need for an 
additional sink and the ground state. 

Relevant environmental parameters going into Lindblad equation ^ are dephasing 
strength 7, recombination rate T and the sink rate n. We use standard values inferred 
from experiments and used before [12], sink rate k = lps' 1 , recombination rate 
T = Ins -1 and dephasing rate at room temperature 7 = 300 cm' ^ Dephasing 
rate, being a product of temperature and the derivative of the spectral density, can 
be estimated by using experimentally determined parameters of the spectral density (as 
done e.g. in reference [12])- This value approximately agrees with the optimal dephasing 
rate at which transfer is most efficient [121 E] • We note that the results shown depend 
very weakly on the actual value of the dephasing rate as long as it is of the same order 



of magnitude as 7 = 300 cm 1 . In Appendix A we show that the values of 7 = 150 cm 



and 600 cm 1 give almost the same optimal PPC size. 
2.2. Overlapping discs model 

Because we want to study the dependence of efficency on the size of PPC, keeping 
the number of pigments fixed, we have to account for the size-dependence of the 
Hamiltonian. On-site energies and interaction strengths in equation ([!]) will be 
determined from the geometry of pigments configuration. Changing the PPC's size two 
gross effects are at play. First, as the distance between pigments is reduced, the dipole- 
dipole interaction between pigments gets larger, enhancing the transfer of excitation 
among them; this effect is already taken into account by the ~ 1/r 3 dependence of 
V nm in equation §2§. Secondly, as chromophores get even closer together, approaching 
the distances comparable to the extent of the chromophore electronic orbitals, the 
Hamiltonian ([!]) is not sufficient for the description of EET anymore because the effects 
due to electronic orbital overlap become important. A detailed analysis of processes 
that take place as chromophores get close together would require advanced quantum 
chemistry methods and is out of scope of this paper. However, the main effect can 
be effectively taken into account by appropriately rescaling parameters of equation 
(JTJ. Because the pigment molecules will be deformed, their excitation energies will 
also change. As the on-site energies e n , being of the order of few eV, are about 
~ 100 times larger than V nm , even a small relative change in e n can have a large 
effect. Effectively, close or even overlapping chromophores will therefore result in widely 
different values of on-site energies e n at different sites, i.e., in a disorder. Thus there 
are two competing factors that determine the optimal size of PPCs: reducing inter- 

|| The values entering the Lindblad equation ^ should be in units of frequency, e.g. s _1 . The 
conversion from inverse centimeters cm' 1 as used traditionally in spectroscopy is given by wis' 1 ] = 
2007rci>, where c = 3x 10 8 and v is value in [cm' 1 ]. 




Figure 1. Illustration of the simple model of PPC. (a) BChl molecule with omitted 
phytyl tail, enclosed by the cylinder as used in simulations (a — 1 A, r = 4 A), (b) 
Discs, randomly positioned in sphere, prior to the reduction of radii of overlapping 
discs. Green discs at the top / bottom represent the input site / output site, (c) Same 
sample as in (b) after the shrinking of overlapping discs. 



chromophore distance increases EET, while the overlapping of chromophores introduces 
a disorder that effectively suppresses EET. 

On-site energy has an additional contribution due to local environment (e.g. 
because of pigment-protein interactions), which is however of the order of ~ 10 2 cm _1 
and is usually much smaller that the on-site disorder due to pigment overlap, which is 
proportional to the unperturbed on-site energy of ~ 10 4 cm _1 . Therefore, the effects of 
pigment-protein interaction were neglected when obtaining the results presented in the 
main text. To verify whether neglecting of on-site disorder due to protein interactions is 
justifiable, we have also calculated the optimal size for iV = 7 pigments with random on- 



site disorder added to each random sample of chromophores. The results (see Appendix 



|B| show that disorder of such magnitudes indeed has no gross effect on the results. 

Each chromophore in our model is represented as a disc - a thin cylinder - of radius 
r and height a. Each disc is supposed to represent an approximate size of the electronic 
cloud of the orbitals involved in the EET (highest occupied electronic orbital, lowest 
unoccupied electronic orbital). We have estimated the height to be a = lA, while 
the cylinder radius is taken as r = 4 A. Size of this cylinder in comparison to a BChl 
pigment molecule can be seen in figure [l^,. Radius r = 4 A is chosen so that it contains 
16 closest non-hydrogen atoms to the Mg atom located in the center of the pyridine 
ring. For given locations and orientations of discs, we then determine if there are any 
overlaps between discs. If two discs overlap, for an example see figure [T]d, we rescale 
the radius of one of them to the new radius r n so that they do not overlap anymore but 
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instead only touch. After eliminating all overlaps we end up with disc's radii r n , for an 
example see figure 

Provided the radius of the n-th disc r n is different from the non-overlapping size 
r, we have to appropriately rescale the on-site energy e n . If the effective size of the 
electronic cloud is reduced from r to r n , the kinetic energy of electron increases by a 
factor r 2 jr\. We therefore estimate that the energy of excitation on a resized disc will 
also scale quadratically with its size, giving the on-site energy dependence 



r(0) 



Tr, 



where is the excitation energy of non-deformed pigments, i.e., the energy difference 
of two lowest electronic states on a pigment. In FMO is approximately ~ 
12 300 cm" 1 . The overall offset of on-site energies is irrelevant for the dynamics in the 
model, therefore we shift all energies by e^. Such quadratic on-site energy scaling 
can be rigorously shown under an assumption that the electronic eigenfunctions of the 
rescaled pigment are just the rescaled eigenfunctions of the original pigment of radius r. 
Let orbitals ipi be the eigenfunctions of the Hamiltonian H(x) = T(x) + U(x), where 
T(x) is kinetic energy operator and U(x) is a confining potential. The on-site energy of 
a given chromophore e n is the difference between the energy of ground and 1st excited 
state, e n = E 2 — E%. Assuming that eigenstates ipi are just scaled to a smaller volume, 
ip*(x) = \ 3 / 2 ipi(x/\), the scaling of on-site energies from equation Q is obtained by 
comparison of eigenvalue equations for the original eigenstate Hipi = E^i and the 
scaled eigenstate H*ip* = E*i/j*, where the scaled confining potential in H* has to be 
U*{x) = U(x/X)/X 2 . 

Dipole strength of the chromophores is similarly scaled linearly with the radius r n 
of the cylinder, 

d n = r -d, (8) 

r 

where d is the bare transition dipole moment of the original chromophore of size r, 
and d n is the scaled dipole strength of the resized disc. This can be justified on the 
same grounds as the scaling of on-site energies, by inserting the rescaled wavefunction 
ip*(x) into the expression for transition dipole matrix of relevant chromophore transition, 
d = (ip g \ex\ip e ). 

There are different possibilities of how to precisely resize discs in order to avoid 
overlaps. While different procedures lead to different on-site energies, the determined 
optimal complex size changes by little. Results in the main text were obtained 
by sequentially inspecting each pair of discs, resizing only the disc having greater 
radius afterwards, while keeping the other disc intact. We have verified other resizing 
procedures, for instance, resizing both discs in pair to the same size. Such resizing 
effectively reduces disorder of on-site energies as even strongly overlapping pigments will 
have identical on-site energies. Nevertheless, for such resizing procedure, the determined 
optimal radius of PPCs are within 2 A of the values obtained by the resizing procedure 
used throughout the paper, and are thus within the error bounds of the model. 
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To summarize, in our overlapping disc model the matrix elements of H are 
calculated for given PPC configuration (positions, as well as disc and dipole orientations) 
by first resizing all overlapping discs, obtaining rescaled radii r n and then scaling dipole 
strengths and on-site energies according to equations ^ and (|8]). 

2.3. Optimality criteria 

We have already introduced equations of motion that govern the dynamics of excitation 
on chromophores, as well the overlapping disc model that allows us to determine the 
Hamiltonian for a given configuration of chromophores. What is left are criteria that 
will enable us to determine whether a given configuration of chromophores is efficient 
in terms of EET. The efficiency of the PPC complex is characterized by the probability 
that the excitation, initially localized on the input site, will be funneled to the reaction 
center trough the output site. For an example of time evolution see |Appendix C The 



efficiency in the model is not unity because the excitation can be lost. The probability 
that the excitation will be transported to the reaction center can be expressed as 

rj = 2K dtp NN (t), (9) 



o 

which will be used as our main efficiency criterion. Closely related is the average transfer 
time, which signifies the speed of transfer of excitation to the reaction center, and is 
expressed as 

2k r°° 

T=— dtt PNN {t), (10) 



with smaller transfer times being better. 

As an additional viability criterion of PPC, robustness of efficiency to static disorder 
will be also inspected. Dynamic disorder due to thermal motion is already effectively 
described by the dephasing Lindblad terms in equation Static disorder due to 

structural changes of PPC, for instance due to changes in biological environment, like 
temperature, electric charges, etc., should be treated separately. A given configuration 
of pigments in PPC is robust to the static disorder if random displacements of pigments 
from the original locations do not induce large changes in PPC's efficiency r\ (or 
equivalently, the average transfer time r). To put it on a more quantitative ground, we 
define the pigment configuration robustness a v (x) for a given configuration of pigment^] 
with positions x = (a?i, a?2, £Cjv)j as a standard deviation of efficiency r\ when pigment 
coordinates are varied in the neighborhood of original positions, 

a ri( x ) = J(V 2 { X + V)~ v{xf)w(y) dy, fj(x) = J i](x + y)w{y) dy, (11) 

Probability density w(y) defines the neighborhood of a given configuration, and is 
localized around the original location of the pigments. The most straight-forward 

% We omitted the disc and dipole orientations from the definition of pigment configuration robustness 
to simplify the expressions. However, no qualitative differences are to be expected if orientations are 
also varied when inspecting the robustness. 
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R [A] R [A] 

Figure 2. Probability density pr,n(v) f° r N — 7 for a range of sphere radii R (density 
plot in the background of the left plot; contours connecting equal values of pr.n are 
also shown). Solid blue curve is the average efficiency (r)(R)). On the right a close-up 
of (77) and the average transfer time (r) (dotted curve, right axis) around the maximum 
of (77) is plotted. 



choice for the distribution w(y) is a product of uncorrelated normal distributions at 
each pigment location, 

w{y) = {2^)- Nl2 ^{ y 0), (12) 

where a defines the size of neighborhood in which the robustness is being probed. 
With given probability distribution, the robustness a v is a function of original pigment 
locations x and size of deviations from original locations a. In the limit of small 
deviations, a — > 0, the expression can be simplified to 

dx'. 



"S(aO = <r 2 $>i(xi) 2 , ruixi) = — (13) 



_1 <-^j ~i ~ Xl 

where the sum goes over all components of pigment coordinates. The robustness a v 
in the limit of small pigment displacements is thus proportional to the amplitude of 
pigment displacements. 
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Figure 3. (a) The optimal enclosing radius R op t(N) of the sphere for different number 
of chromophores N, based on the maximum of the average efficiency (r)(R)). The width 
of the shading (i.e. « ±lA) denotes a range of R for which the average efficiency is 
within 1% of maximal value. Vertical dashed line marks the case of FMO with N = 7 
chromophores. (b) Structure of FMO complex as determined from spectroscopic data, 
enclosed into the optimal sphere of the radius R ~ 16 A. The structural data was 
obtained from PDB entry 3EOJ [3] 

3. Optimal number of pigments: the case of Fenna-Matthews-Olson 
complex 

In this section, the efficiency and robustness of random configurations within the model 
will be considered on an example of FMO complex. FMO consists of N = 7 BChl 
pigments. On-site energy for BChl was chosen to be = 12 300 cm -1 , which is within 
the range of on-site energies for BChls in FMO as determined in the literature (221E31I21]- 
The strength of transition dipole moment d — \d\ was taken as d 2 = 26 D 2 [^] (note that 
published values for d from calculations and experimental data vary considerably [315]). 
On-site energy and dipole strength d used hold for BChl pigments in general and 
therefore the results presented are expected to be valid also for other PPCs containing 
BChls, not just for the FMO complex. 

To determine the optimal size of PPC for a given number of chromophores (or 
equivalently, optimal number of chromophores for a given size), we considered two 
criteria based on overall behavior of efficiencies and robustness of random configurations. 
In the following, the motivation for choosing such optimality criteria will be given, and 
the results for the case of FMO will be presented. 

+ Dipole moment is given in units of debye (D). Conversion of interaction V mn to units of cm -1 is 
obtained by conversion jEL_ — 5030 cm -1 A 3 . 
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3.1. Average efficiency 

We shall use the average efficiency (rj) , averaged over random positions and orientations 
of pigments enclosed in a predefined volume. The reason to use random averaging 
with uniform distribution is twofold: first, high average EET efficiency under 
uniform averaging will mean that there are many different configurations that have 
high efficiency, i.e., high efficiency is globally robust. Choosing averaging over 
random configurations therefore offers insight in how special efficient configurations 
of chromophores are within the space of all configurations. Second reason is that we 
have a priori no knowledge what would be the appropriate measure for possible pigment 
configurations under say different protein cage foldings due to for instance mutations. 
A uniform measure represents in this case a "least-information" distribution. F] Using 
configurations sampled according to a uniform distribution over chromophore positions 
within a ball of radius R and random orientations of dipoles and discs, the average 
efficiency is calculated. Formally, it can be written as 



where X contains positions and orientations of chromophore discs and dipoles (apart 
from the positions of input and output sites which are fixed on the poles of the sphere), 
and u> con f(X) oc 1 signifies a uniform distribution of chromophores inside the sphere. 
Observing the dependence of the average efficiency (t](R)) for different number of 
chromophores and different radii R of the enclosing sphere, we can determine the optimal 
number of chromophores for a given radius R, or equivalently, the optimal radius R opt 
for a given number of chromophores. 

To obtain a more detailed information about efficiencies of random configuration, 
we also observed the probability distribution over efficiencies Pn t N(ji)i defined as 
Pr,n{v) = Ir^(v(^) ~ v) w conf(X)dX . For the number of chromophores as found in 
FMO (N = 7), the probability distribution over efficiencies Prj{j]) is shown in figure [2] 
When going from large radii R to smaller, configurations tend to get more efficient, 
which is expected as the chromophores are closer to each other, thus increasing the 
dipole coupling. However, as R is reduced even further, overlapping of chromophores 
gets more probable, causing an on-site energy disorder. This leads to the localization 
of excitation on chromophores not connected to the sink site. Such configurations have 
low efficiency. Therefore, as R gets smaller the distribution Pr^n becomes bimodal, with 
lower efficiency mode due to overlapping configurations and high efficiency mode for 
non-overlapping configurations. 

Low efficiency of overlapping configurations therefore leads to a maximum in the 
average transfer efficiency (i](R)) at the optimal radius -R pt- The average transfer 
efficiencies at the optimal radius are rather high, e.g. for N = 7 in figure [2] it is 

* We have to note that we also checked other distributions, for instance a uniform distribution on the 
surface of a sphere of radius i?, and obtained practically the same results. For instance, the difference 
in the position of the maximum in figure [2] was within our error estimate of 1 A (seen as an "error" 
band in figure |3|. 




(14) 
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(i](R opt )) pa 0.95 with a large fraction of configurations having even larger efficiency 
than the average. Thus within the model, high efficiency is not due to finely tuned 
pigment positions and orientations, but occurs for majority of pigment configurations 
for parameters estimated to be relevant in PPCs. For the FMO case with N = 7, 
the optimal radius was estimated to R opt ~ 16 A, which fits the actual configuration 
of pigments very well (see figure [3]). The average transfer time (r) is also minimal at 
R = R opt (see figure [2]). Optimal average transfer time of ~ 30 ps is so large due to the 
contribution of very inefficient configurations of chromophores. Looking at the average 
transfer time of the 5% of most efficient configurations, we get the value of 5 ps, which 
is comparable to the transfer times as determined using different models of the FMO in 
references [El EH El ES] ■ 

The estimated optimal radius is quite insensitive to small variations of input 
parameters, e.g. dipole moment d, chromophore disc radius r or its thickness a, or the 
scaling of on-site energies and dipole strengths of resized discs. For instance, decreasing 
disc radius to r = 3.5 A decreases R op t by ~ 2 A, changing disc thickness to a = 0.5 A 
or 1.5 A changes i? pt by ~ =F2 A, while changing quadratic energy scaling to a linear or 
cubic one again changes R opt by « =f2 A. Similarly, changing the dephasing rate 7 by a 



factor of 2 changes _R pt by ~ 2 A, see Appendix A Details of the disc resizing procedure 



also change R opt for less than 2 A, as the extreme case of resizing each overlapping disc 
pair to the same size reduces R opt by ~ —2 A. 

The optimal radius of the enclosing sphere was obtained from the average 
efficiency over all random configurations within a sphere. However, even though the 
evolutionary drive to more efficient configurations might not be very strong if majority 
of configurations are already efficient, still some optimization is to be expected. Thus 
one might argue that the optimal enclosing volume of the natural PPCs should be 
determined considering only the ensemble of more optimal configurations. We will 
denote such averages with (i]) p where p specifies a portion of most efficient configurations 
that should be taken into account when calculating the average (e.g. (77)0.05 is the average 
of 77 over 5% of most efficient configurations as shown in figure [1]). As p is reduced, 
the overall value of average efficiency (r/) p will increase. The increase will be more 
pronounced in the region of R < R op t, where the distribution is bimodal. The location 
of the maximum of (r/) p will be thus moved to smaller values of R, indicating more 
densely packed chromophore configurations. However, as we will see in next subsection, 
robustness of such densely packed configurations deteriorates very quickly, supporting 
our choice of estimator for the -R op t- 

Note that the overlaps between pigments and protein cage are not considered 
explicitly in the model. If overlaps with protein cage would be taken into account, 
-Ropt would represent the size of a protein cage, whereas in our model without pigment- 
protein overlaps, _R op t is the size of a sphere that contains all pigment centers. For 
instance, looking at figure^, we can see that the sphere with R opt contains all pigment 
centers, while parts of few pigments still protrude the bounding sphere. If overlaps of 
pigments with the protein cage would be taken into account explicitly, i? op t would be 
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Figure 4. The average robustness of top 5% of most efficient configurations (u^o.os 
(solid blue curve) as a function of the enclosing sphere radius R. Dash-dotted green 
curve represents the average robustness of the top 15 % of most efficient configurations 
(tr I? )o.i5. Dotted red curve is the average efficiency over the top 5% of configurations 
(v)o.05 (right axis). 



approximately by a disc radius r = 4 A larger, i.e. in corresponding figure, the bounding 
sphere would enclose all pigments completely. 

3.2. Robustness 

Robustness of PPC configurations to static disorder should also be taken into account 
when determining whether a given configuration of pigments is feasible, as the conditions 
in which PPCs operate are subject to constant environmental changes. In previous 
subsection, we have inspected the probability distribution of efficiencies rj over random 
configurations, showing that majority of random configurations achieve relatively high 
efficiency when the enclosing volume is optimal. In this subsection, we will present 
analogous analysis of the robustness of random configurations, in particular of those 
with high rj. We shall show that highly efficient configurations in small enclosing R are 
very fragile. 



We have defined robustness of efficiency o r] in equation (11). In simulations we 
have displaced the pigment positions according to normal distribution with a width of 
a = 0.1 A, which is small enough to quantify the robustness in the neighborhood of 
specific configuration, while larger than displacements due to thermal vibrations, which 
are already effectively described by the Lindblad equation. We are specifically focusing 
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on a subset of the most optimal configurations in terms of rj. The average robustness 
of the subset of optimal configuration is denoted as (cr^) p , where p specifies a fraction 
of most optimal configurations in terms of EET efficiency 77. As an example, we will 
consider robustness of top 5% of efficient configurations (cr^o.os- The dependence of the 
average robustness on the radius of the enclosing sphere R is shown in figure |4j The 
average efficiency of optimal configurations (77)0.05 is also shown in the figure. While the 
average efficiency of top 5% of optimal configurations continues to rise as the enclosing 
sphere radius R is reduced, we can see that the average robustness (a v ) worsens very 
quickly as the R drops below -R opt . 

Quick worsening of EET robustness with reducing sphere radius suggests that 
even if the PPC configurations occurring in nature are indeed optimized in terms of 
pigment positions and orientations, the excessive stacking of pigments is not favored 
as it makes PPC configurations very sensitive to any displacements of pigments. The 
transition from robust to non-robust regime takes place at a radius comparable to R opt 
at which the average efficiency (7/) has a maximum. This is not surprising as both, the 
efficiency of configurations and robustness of configurations, are strongly influenced by 
the overlapping of pigments which gets more pronounced for R < -R pt- 

We have presented results for the robustness a v of the 5% of most efficient 
configurations, with pigment displacements a — 0.1 A. General characteristics of (cr r? ) p 
however do not quantitatively change for different p (the case of p = 0.15 is also shown in 
figure [4]) or displacements a. Most importantly, the radius R at which the robustness of 
configurations drops significantly takes place at approximately the value of i? pt- Same 
behavior of robustness is observed also in the limit of infinitesimal robustness from 
equation (13) where a — > 0. 



4. Optimal pigment numbers in other PPCs 

In previous section we have calculated the optimal size R or the optimal number 
of pigments for the FMO complex. We also demonstrated that a large portion 
of chromophore configurations has high efficiency when the enclosing volume is 
properly chosen (~ -R pt)- Additionally, robustness of configurations to chromophore 
displacements starts to deteriorate quickly once the enclosing volume is reduced below 
-R opt . Based on these two observations, we argue that the enclosing volume of PPCs 
occurring in nature should be close to the optimal volume as determined by our simple 
model. In this section we will present similar results for the PPCs containing chlorophyll 
(Chi) chromophores. 

We compare results of the model to the structure of PPCs from the Photosystem II 
(PSII) [3Zj, found in cyanobacteria, algae and plants. PSII consists of multiple functional 
units, which are either part of the outer light-harvesting antenna or the inner core, 
to which excitations are funneled. In the light-harvesting antenna we will consider 
the light harvesting complex II (LHCII), while in the core we will focus on the PC43 
and PC47 complexes that funnel excitations to the reaction center and thus have a 
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similar role as the FMO complex in bacteria. A monomeric unit of LHCII contains 14 
chlorophyll molecules (8 Chl-a and 6 Chl-6), while CP43 and CP47 contain 13 and 16 
Chls respectively. 

Model parameters for the sink rate, dephasing and recombination are kept the 
same as in the FMO case, while the transition dipole strength and on-site energies 
are different for Chi molecules. Transition dipole moment of Chi molecules is chosen 
as d 2 = 15 D 2 and the on-site energy e( ) = 15 300 cm -1 , where values were taken 
according to reference [HE] (we take the average between values for Chl-a and Chl- 
b). For the CP43 and CP47 complexes we have simulated random configurations of 
13 and 16 chromophores enclosed into a sphere as the actual chromophore positions 
are distributed relatively uniformly in all directions. The shape of LHCII is however 
significantly elongated in one direction. We therefore choose the cylindrical volume, 
having only one additional parameter that has to be provided, i.e. the ratio between the 
cylinder radius R c and cylinder height A. Based on positions of the LHCII chromophores 
we have estimated the ratio of the two to be R c /A = 0.34. 

The CP43 and CP47 primarily play a role of an exciton wire, making the model 
with input site and output site at the opposite sides of the sphere applicable. The 
optimal radius as predicted by the model is R ~ 18 A for CP43 and R ~ 20 A for CP47. 
As the LHCII also has to transport excitations from adjacent complexes, we have also 
determined the optimal shape of LHCII with input and output sites located at the 
opposite sides of the enclosing cylinder. With the ratio R c /A fixed, we have varied the 
height A of the cylinder and determined the optimal height to be A opt ~ 43 A. 

In addition to acting as an excitation wire, CP43 and CP47 complexes are also 
directly involved in the absorption of photons, in which case the role of the input site 
can be taken by any chromophore site. This is even more common scenario in the LHCII 
complex, whose primary role is the absorption of photons. To verify whether the findings 
about optimal enclosing volume are also valid when the main purpose of PPC is the 
absorption of photons, we randomly placed the input site inside the enclosing geometry 
for each configuration in random ensemble. General characteristics of the distribution 
over efficiencies PB,{rj) do not change considerably, however, the distribution is somewhat 
shifted to the higher efficiencies because in many random configurations the input site 
is considerably closer to the output site than the diameter of the enclosing volume. This 
results in the optimal size of enclosing volume being somewhat larger, R opt 20A for 
CP43 and R opt ~ 22A for CP47. For the LHCII we have moved the output site to the 
midpoint on the side between top and bottom of the cylinder, where the actual output 
site is supposedly located jUJ. For such geometry and previously used R c /A = 0.34, we 
have obtained the optimal height of the enclosing cylinder at A rs 47A. 

The optimal enclosing values obtained from the model (averaged between the case 
for fixed input site and random input site) were compared to the actual configurations 
of pigments as obtained from spectroscopic data, and are shown in figure |4^,-c, showing 
good agreement. For the spherical geometries, we have centered the sphere of optimal 
radius R op t to the arithmetic mean of locations of BChl/Chl centers. For the LHCII, 
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(a) CP43, R opt = 19 A 



(b) CP47, i? opt = 21 A 





(c) LHCII, A opt = 45 A 



(d) 
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Figure 5. (a) - (c) Various pigment-protein complexes (PPCs) enclosed in optimal 
geometries as described in the main text. The structural data was obtained from PDB 
entries 1RWT (LHCII) 39J and 3 ARC (CP43, CP47) 40J. In (c) enclosing geometries 
of two additional monomeric units of the LHCII trimer are also shown, (d) Plot 
showing dependence of the average efficiency (?y) on the size of the enclosing geometry. 
The upper panel shows the case of CP43 (green) and CP47 (blue), while the lower panel 
shows the case of LHCII. Solid curves are for the fixed input site, while dotted line 
for the randomly placed input site. Vertical dashed lines mark the sizes of enclosing 
volumes as used in subfigures (a) - (c) . 



where cylinder was used as the enclosing geometry, the cylinder axis was determined 
such that r ±i wa s minimal, where r±i is the distance from the cylinder axis to the 
position of i-th Chi center. Interestingly, three cylinder axes do not lie in a plane but 
are instead tilted at an angle 15° to the plane containing three cylinder centers. It is 
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not known if this plays any functional role. 
5. Conclusion 

We have studied the efficiency of excitation energy transfer in protein-pigment complexes 
for random configurations of pigments. The Hamiltonian part of Lindblad master 
equation is determined from the geometry of pigment configurations. If pigments are 
too close, so that they overlap, this introduces a disorder in on-site energies, effectively 
inhibiting excitation transport. Fixing the enclosing volume in which pigment molecules 
are located we have calculated the average efficiency over random pigment configurations 
as well as robustness of efficiency to variations of pigment locations. Doing this we have 
determined the optimal number of pigments for a given size of the complex. Even 
though the model is an oversimplification of actual processes that take place in nature, 
statistical predictions obtained from the model are robust to its variations. 

Comparing theoretically predicted optimal number of pigments with several 
naturally-occurring complexes we find good agreement. This might indicate that 
PPCs are not optimized just to have the highest possible efficiency - in fact, efficient 
configurations are quite common - but instead to be robust to variations in pigment 
locations. Namely, it turns out that configurations optimized for the highest efficiency, 
that is those with specially tuned positions and dipole orientations, are very sensitive 
to small perturbations. The number of pigments in nature is therefore chosen in such 
a way that the probability of having efficient configurations that are at the same time 
also robust is the highest. 

The presented findings could be in principle verified experimentally by modifying 
the structure of known PPCs and probing the efficiency of excitation transfer. For 
the FMO complex the structure was already changed by mutation of genes encoding 
the structure of BChls, as well as by substituting the carbon 12 C atoms with 13 C [22J. 
Comparison of excitonic spectra revealed no distinctive differences in the dynamics of 
excitations, complying with the hypothesis that configurations are not highly tuned but 
instead very robust. An additional intriguing possibility would be also to inspect the 
characteristics of FMO with mutated protein cage, modifying positions and orientations 
of pigment molecules. One could also compare our predictions for the optimal sizes (e.g. 
figure [3k) with other complexes occurring in nature. 




Appendix A. Dependence on the dephasing rate 7 

In the simplified model used, environmental interaction is described by dephasing rate 7, 
having the same value for all chromophore sites. In principle, environmental interaction 
requires more involved equations of motion (e.g. HEOM [7J), taking into account the 
spectral density of the environmental modes. However, due to crude nature of the model, 
simplified description is expected to account for main environmental effects influencing 
the efficiency of exitation transfer (see e.g. [501 I3"T] for more detailed comparison 
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Figure Al. (a) Average efficiencies (77) for different values of dephasing 7 = 
150 cm -1 , 300 cm" 1 , 600 cm -1 , (b) Average efficiencies (rf) for different values of 
random on-site energies disorder er ra nd- 



of approaches). The adequacy of simple Lindblad-type description of dynamics for 
determination of optimal size is also justified due to the high robustness of the results 
to the actual choice of dephasing value 7, as seen in figure ATi, where R opt only changes 
for ~ ±2A as dephasing rate 7 is changed by a factor of 2. Optimal size R opt of the 
PPC is somewhat smaller as dephasing rate gets stronger, which is expected as larger 
dephasing rate enables transfer across the sites with greater on-site energy mismatch, 
getting increasingly common in more compact configurations of chromophores. 



Appendix B. Random on-site disorder 

To verify whether effects of the local chromophore environment due to e.g. pigment- 
protein interaction can affect the findings about the optimal PPCs sizes, we have 
amended the Hamiltonian in equation ([TJ with random on-site disorder e™ nd of the 
magnitudes as found in naturally occurring PPCs (i.e. on-site energy differences in the 
order of ~ 100 cm -1 ). The values of disorder for each realization of random PPC were 
calculated according to Gaussian distribution with variance <7 ra nd- Results are shown in 
figure AT d. In the region R < R opt , where average transfer efficiency is strongly affected 



by disc overlaps, an addition of random on-site energy disorder has no noticeable effect. 
The effect is more pronounced for R > R opt where overlapping of discs is not the limiting 
factor of transfer efficiency any more. The estimated optimal size i? op t however is not 
changed considerably by an addition of random on-site disorder. 



Appendix C. Time evolution of site populations 

To provide some insight into the temporal dynamics of the excitation transport, we 
present the time evolution of site populations for two different realizations of PPC 
within the Lindblad model. In figure [CTJ t, time evolution for the FMO Hamiltonian 
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Figure CI. (a) Time evolution of site populations (n\p\n), calculated using Lindblad 
model for the FMO Hamiltonian from reference [5], resulting in efficiency of r\ = 0.99 
and average transfer time of r = 6.2 ps. At t — 0, excitation is localized on site 1. Sink 
is connected to the site 3. Dotted line represents the population of the sink, (b) Time 
evolution for a random sample, generated for R — 18 A, with efficiency rj = 0.90 and 
average transfer time r = 45 ps. Blue line is the population of the input site that is 
initially excited. Red line is the population of the output site, connected to the sink. 
Dotted line represents the population of the sink. 



from reference [B] is shown, and in figure CI d, the time evolution of randomly generated 
PPC of radius R = 18 A. The values of dephasing, sink rate and recombination rate are 
the same as used in the main text. 
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